qbit=1;showfig = 0;

Nsub=256;   % no. of subcarriers for rate computation
Ngrid = 1000;
L = 10e-1;

% factor = 1.6;
% x1 = 40e-1*factor;
% y1 = 50e-1*ones(size(factor));
% z1 = 90e-1*factor;
% x2 = -30e-1*factor;
% y2 = 50e-1*ones(size(factor));
% z2 = 40e-1*factor;
% TRX = zeros(6,length(factor));
% TRX(1,:) = x1;
% TRX(2,:) = y1;
% TRX(3,:) = z1;
% TRX(4,:) = x2;
% TRX(5,:) = y2;
% TRX(6,:) = z2;
load('TRX_std.mat','TRX','NTRX')
NTRX = 200;
dist = rms(TRX(1:3,:))*sqrt(3)+rms(TRX(4:6,:))*sqrt(3);

fc = 30e9;
c = 3e8;
BW_des = 1.5e9;
Ptx = 10;%dBm
nBS = round(2.^(0:0.5:5));
nBS = unique(nBS);
% nBS = [1 2];

cal_rate = 0;
if cal_rate

Rates = zeros(4,length(Ptx),NTRX);
for n = 1:length(nBS)
    n1 = nBS(n);
    n2 = nBS(n);

    parfor m = 1:NTRX
        BS = struct();
        BS.n1 = n1;
        BS.n2 = n2;
    bsxyz = TRX(:,m);
    bsxyz = bsxyz(1:3);
    D = [bsxyz randn(3,2)];
    [Q,~] = qr(D);
    BS.d1 = Q(:,2)';
    BS.d2 = Q(:,3)';
   Rates(:,n,m) = MISO_compute_rate(10,Nsub,L,TRX(:,m),fc,BW_des,Ptx,BS,showfig); 
    % Rates1(:,n,m) = my_compute_rate(1,Nsub,L,TRX(:,m),fc,BW_des,ptx,showfig)/BW_des; 
    % Rates2(:,n,m) = my_compute_rate(2,Nsub,L,TRX(:,m),fc,BW_des,ptx,showfig)/BW_des; 
    % Rates3(:,n,m) = my_compute_rate(3,Nsub,L,TRX(:,m),fc,BW_des,ptx,showfig)/BW_des; 
    end
    disp(n)
end
Rates = mean(Rates,3);
save('data/Rates_sum_rate_vs_nBS.mat','Rates','Ptx')
else
    load('data/Rates_sum_rate_vs_nBS.mat');
end

set(0, 'DefaultAxesFontName', 'Times New Roman');
semilogx(nBS.^2,Rates(2,:),'bs-','LineWidth',2.5,'MarkerSize',10) % optimal
hold on
% semilogx(nBS.^2,Rates(2,:),'v-','LineWidth',2.5) % proposed
semilogx(nBS.^2,Rates(3,:),'r^-.','LineWidth',2.5,'MarkerSize',10) % proposed TTD 
semilogx(nBS.^2,Rates(4,:),'k+--','LineWidth',2.5,'MarkerSize',10) % proposed TTD 
xticks(nBS.^2)
h=plot(dist(round(end/2)),0);
h = h.Parent;
xlabel('Number of basestation antennas \it{N}^{BS}','FontSize',16)
ylabel('Average achievable rate (bit/s)','FontSize',16)
legend('FZ-GSA, actual BS','FZ-GSA, approximated point source' ...
    ,'Optimal beamforming' ...
    ,'Location','southeast','FontSize',16)
hold off
set(get(gca, 'Legend'), 'Box', 'off'); 
set(gca, 'FontSize', 16);
set(gca,'FontName','Times New Roman')
set(gcf, 'Position', [100 100 600 500]); % 调整当前图窗的大小

saveas(gcf,'fig/Fig_sum_rate_vs_nBS.fig')
saveas(gcf,'fig/Fig_sum_rate_vs_nBS.png')